#------------------------------------------------------------------------------
# Import libraries
#------------------------------------------------------------------------------

rm(list=ls())
library(LalRUtils)
libreq(tidyverse, data.table, zoo, tictoc, 
       fst, fixest, PanelMatch, patchwork,
       rio, magrittr, janitor, did, panelView, 
       ggiplot, tictoc, binsreg, interflex)
set.seed(42)
theme_set(lal_plot_theme())

#------------------------------------------------------------------------------





#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])

#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
setnames( vcf, "d", "D")
vcf_data <- copy( vcf )
gfc <- fread("gfc_dta.csv" )


#------------------------------------------------------------------------------



#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

# %%
vcf = vcf_data[year>= 1990]

# %% one-shot treatment in 2008
vcf[, pref_bin := ntile(cover_1990, 10)]
vcf[, time_fra := year - 2008]
vcf[, D_fra := sch * (year >= 2008)]
fitter = function(cut) {
  m = feols(forest_index ~ D_fra | cellid + styear, 
            data = vcf[pref_bin >= cut], 
            cluster = ~blk)
  tidy(m)[, 2:3]
}

cutmods = map_dfr(1:10, fitter) %>% setDT
cutmods$n = 1:10
colnames(cutmods)[1:2] = c('beta', 'se')


# %% GFC estimate for PRI
gfc[, D_fra := sch * (year >= 2008)]
fitter = function(cutoff){
  dat = gfc[pref_bin >= cutoff & gfc3 == TRUE]
  m = feols(def_ha ~ D_fra | village[t] + styear, cluster = ~block, dat)
  tidy(m)[, 2:3]
}

tic()
cutoff_res = map_dfr(1:10, fitter) %>% setDT
toc()
cutoff_res$n = 1:10
colnames(cutoff_res)[1:2] = c('beta', 'se')
# %%
(FRA_fit_gfc = ggplot(cutoff_res, aes(n, beta)) +
    geom_point(size = 2) +
    geom_errorbar(aes(ymax = beta + 1.96 * se,
                      ymin = beta - 1.96 * se), 
                  alpha = 1, width = 0) +
    geom_hline(yintercept = 0) +
    scale_x_continuous(breaks = 1:10) +
    labs(title = 'GFC', y = "Effect (Deforested Aarea)",
         x = "Inclusion Threshold Decile of Forest Cover in 2000 ",
         caption = '') + 
    theme( legend.position = 'none', 
           axis.title = element_text( size = 23 ), 
           axis.text = element_text(size = 23), 
           plot.subtitle = element_text(size = 25))
)
# %%
ggsave( "appendix_figureA3.pdf", height = 5, width = 12, 
        device = cairo_pdf)
# %%

#------------------------------------------------------------------------------




